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Abstract 

The simulation of quantum systems is a task for which quantum com- 
puters are beheved to give an exponential speedup as compared to 
classical ones. While ground states of one-dimensional systems can 
be efficiently approximated using Matrix Product States (MPS), their 
time evolution can encode quantum computations, so that simulating 
the latter should be hard classically. However, one might believe that 
for systems with high enough symmetry, and thus insufficient param- 
eters to encode a quantum computation, efficient classical simulation 
is possible. We discuss supporting evidence to the contrary: We pro- 
vide a rigorous proof of the observation that a time independent local 
Hamiltonian can yield a linear increase of the entropy when acting 
on a product state in a translational invariant framework. This cri- 
terion has to be met by any classical simulation method, which in 
particular implies that every global approximation of the evolution 
requires exponential resources for any MPS based method. 



1 Introduction 

One of the main challenges in the field of Quantum Computation is to deter- 
mine the boundaries between classical and quantum computers and pinpoint 
the computational problems for which quantum devices will be more powerful 
than classical ones. The research on this topic goes in two directions: On the 
one hand, several efficient quantum algorithms for mathematical problems 
which are believed to be hard classically have been found, the most prominent 
being Shor's celebrated algorithm for factoring and discrete logarithms [1]. 
On the other hand, the idea to use quantum simulators to simulate quan- 
tum systems has gained growing attention, as simulating quantum systems is 
arguably the most natural thing to do with a quantum computer [2]. More- 
over, in the near future quantum devices are far more likely to be used for 
simulating quantum systems than for factoring large numbers, both due to 
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the much smaller number of qubits required and the weaker demands on the 
control of the system. In particular, there is no need to have a universal 
quantum computer available [3-5]. 

In fact, while universal quantum computation requires a dynamics which 
is either inhomogeneous in time [6, 7] or in space [8, 9] (or both, as in the 
standard network model), physically interesting dynamics to simulate can 
be homogeneous in time, i.e., generated by a time- independent Hamiltonian, 
and translational invariant. 

Regarding static properties translational invariance appears to consid- 
erably simplify the problem and enables us in many cases to find efficient 
classical approximation algorithms. For instance ground states of local one- 
dimensional Hamiltonians, in particular in the presence of an energy gap, 
have an efficient classical description [10, 11] in terms of so-called Matrix 
Product States (MPS) [12], which are used as a variational ansatz in the the 
Density Matrix Renormalization Group (DMRG) algorithm [13,14]. 

If static properties of a quantum system are efficiently simulatable on a 
classical computer in the presence of translational symmetry, why not also 
time evolution governed by a time-independent local Hamiltonian? 

The present paper is devoted to providing evidence for the hardness of 
classically simulating quantum mechanical time evolution for systems which 
are homogeneous in time and space. Hence, quantum computers (or simula- 
tors) may indeed yield an exponential speed-up compared to their classical 
counterparts. 

Clearly, we cannot hope for a complexity theoretic separation between 
classical an quantum due to the lack of parameters needed to encode a com- 
putational task. We will rather show first that a sufficiently fast (e.g. linear) 
increase of the entanglement entropy makes it hard - and for MPS based 
approaches [15, 16] impossible - to even approximate the state of the sys- 
tem efficiently. That such a linear increase can occur even in the simplest 
one-dimensional systems was observed in [17] based on arguments of confor- 
mal field theory together with numerical diagonalization of chains up to 100 
sites. The main part of the present paper is devoted to a rigorous proof of 
this observation. 

2 The hardness of simulating time evolution 

We show the hardness of simulating time evolution using MPS in two steps: 
In a first stage, we take a particular system initially in a translational invari- 
ant product state of spin-| particles, and let it evolve under a translational 
invariant nearest neighbor Hamiltonian. For this system, we prove a linear 
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lower bound to the entropy of any contiguous block as a function of time, 
as observed by Calabrese and Cardy in [17]. This evolution has also been 
studied in [18], where similar effects were observed. We then combine this 
linear bound with the continuity inequality for the von Neumann entropy 
which leads to a linear lower bound on the block entropy of any approxi- 
mation of the state [19]. As the computational effort needed to deal with 
a general MPS grows exponentially in the entropy, it would thus increase 
exponentially with the time t to be simulated. 

In the following, logarithms are understood to the base 2. In particular, 
we define the von Neumann entropy S{p) = — tr[plogp] with the base 2 
logarithm. The natural logarithm will be denoted by In. 

2.1 Setup and entropy scaling 

Theorem 1 Consider a chain with an odd number N of spins with periodic 
boundary conditions, corresponding to a Hilbert space (C^)®^, which at time 
t = is in the state \4'{0)) = | | . . . |) = |1 . . . 1) and evolves under the Ising 
Hamiltoniai^ 

so that \^{t)) = e-''^^\ip{0)) . Define 

SLit) = -tT[pL{t) log PL (t)] 

With pL{t) = tTL+i^,„^N\i;{t)){ij{t)\. Then for N> 20, L > 10, 4 < t < eL/4, 
and t < N/5, 

SL{t)>^t-l\nt-l . 

We postpone the proof of the Theorem to Sec. [3l 

Note that although the lower bound is independent on L - the intuition 
behind being that entanglement arises at the boundaries - the maximally 
attainable entanglement is proportional to the length of the block due to the 
constraint t < eL/A. 

2.2 Exponential scaling of bond dimension 

Let us now see what the entropy scaling implies for an ansatz state |0) = 
\(f>{t)) which we use to approximate \ip) = \ipit)). Let pL = Phit) = tYB\ip){4'\y 

^ We choose (j^, <Jz the Pauh matrices with eigenvalues ±1. 
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ai = CTLii) = tTB\(p){(p\ be the reduced states on a block A of length L {B is 
the complement of A). The error we make in this approximation is 

e>m){^\-\(t>){mi>\\pL-(rL\\i=: 2T , 

where we have used that the trace norm ||?7||i = tr |?7| is contractive under 
the partial trace. We now exploit Fannes' continuity inequality for the von 
Neumann entropy in its improved version by Audenaert [20], 

\SipL) - <Tlog(2^ - 1) + HiT, 1 - T) 

<TL + 1 , 

where H{T, 1 — T) = —T log T — (1 — T) log(l — T) < 1 is the binary entropy, 
and thus 

SicTL) > SipL) -TL-1 
> S{pl) -\eL-l. 

Using the lower bound on S{pi(t)) = Siit) of Theorem [T], we have (under 
the corresponding conditions on A^, L, and t) that 

>^t-leL-^\nt-2. 
on 

To get an optimal bound, we choose L as small as possible for the given t, 
L = 4:t/e, which implies the new constraint t > 5e/2, and obain 

^KW)> (^-|)t-ilnt-2. (2) 

In order to turn this into a lower bound on the size of an MPS 
used to approximate the state at time t, note that the entropy of any block is 
naturally upper bounded by the rank of the reduced state pi, which in turn 
is at most D^. Here, D is the so-called bond dimension of the MPS, which 
is polynomially related to the complexity of the MPS [14]. Then, ([2]) gives 

logZ^> (l^-^)t-ilnt-l, (3) 

i.e. the bond dimension - and thus the resources for a straight forward MPS 
simulation of the time evolution - will scale exponentially with time as soon 
as the accuracy e < eo = 2e/37r ^ 0.577. 

Note that the same argument holds similarly for lattices in higher spatial 
dimensions. As long as the entropy S{pA(t)) of a region A has a linear 
leading term, there is an accuracy eo beyond which S{aA(t)) has to grow 
linearly as well. Hence, approaches based on projected entangled pair states 
(PEPS) [21] or variants thereof are burdened with the same problem. 
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3 Proof of Theorem [T] 



We now prove Theorem [H Let us briefly sketch the main steps: First, we 
map the spin chain to a one-dimensional system of free fermions, which can 
be solved exactly, and bound the error made by going to the thermodymamic 
limit ^ oo. The entropy of a contiguous subblock of length L is equal 
to the entropy of the corresponding fermionic modes; using a parabola as 
a bound on the binary entropy we can lower bound the block entropy by a 
quadratic function in the correlation matrix. Combining this with the exact 
solution for the thermodynamic limit, we obtain a lower bound on the block 
entropy which is essentially of the form 

with Jn Bessel functions of the first kind. In the appendix we bound this 
sum by deriving a lower bound for L = oo and an upper bound on the error 
made by extending the sum. Taking all this together then proves Theorem [H 
Before starting with the proof, let us recall the conditions imposed, namely 
> 20, L > 10, 4 < t < eL/4, and t < N/5. In addition we can assume 
L < N/2 as the overall state is pure. 

3.1 Diagonalization of the Hamiltonian 

We start by mapping the spin operators to fermionic Majorana operators via 
the Jordan- Wigner transformation 

n—l n—l 

Cn = l[cr^^^cTi!^K C+AT = n n = 0,...,N-l. (4) 

1=1 1=1 

We will refer to Cq, . . . , c^^i as the position and to cn, . . . , C2n-i as the 
momentum operators; each pair Ck, CN^k defines a fermionic mode. For the 
initial state \ip{0)) = |1 ■ ■ ■ 1) we have 



AT 



1^(0)) = TT i(l + ickCk+N) = lim 



CkCk+N 



k=l 

and the Hamiltonian 

n = 



fs^oo 2cosh(/3) 



^^(Cj — Cj+1 mod N)Cj+N 



(5) 
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is equal to the spin Hamiltonian ([T]) of Theorem [T] on the —1 eigenspace of 



n = ai^^ ■ ■ -aT' 1 while for the +1 eigenspace, the coupling between spins N 
and 1 has opposite sign. Since n|'?/;(0)) = —\ip{Q)) for odd and \H, 11] = 0, 
([5]) indeed describes the evolution of the Hamiltonian of Theorem [H Note 
that the following results equally hold for even where they describe an 
Ising system with a flipped coupling across the boundary. 

Since the initial state is a Gaussian and the Hamiltonian a quadratic 
function in the Majorana operators, the system remains in a Gaussian state 
throughout its evolution and is thus completely characterized by its second 
moments, the antisymmetric correlation matrix [22,23] 



AN) 



\i = -ftr (p[cfc,Q]) 

0-1 
1 

J- 1 — e i oe , 
where we defined the Hamiltonian matrix H by 

. 2Af-l 



which for the initial state is 



It evolves as 



(6) 



= - ^ Hki[ck,ci] . 



k,l=0 



We now apply Fourier transformations J-'ki = e^'^'^^^l^ j \pN to both the posi- 
tion and the momentum subspace. Thereby, Fq and H become block diagonal 
with blocks 

fo(0.)=f? ~}\ (7) 



and 



m<\>k) 



1 



1 - e*<^^ 



where 0^ = 27r-^. The evolution ^ becomes 



exp 



-Hi 



fo(0fc) exp 



(9) 



such that the CM at time t (written in modewise ordering) is 

/ 7o 7i ■ ■ ■ 77V-1 \ 



7-1 7o 



V 71- 



(10) 



■TV 



70 / 
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with 



Af-l 



N ^ 



fc=0 



Solving ([H]), we find that 



In 



in 9n 
9~n fn 



with 



7V-1 



^" = ^ [e-"^^l^ + e'"^^'^] sin(2t sin f ) 



fc=0 



'IV. 



(12) 



^ [1 - e'"^' + (1 + e'"^') cos(2t sin f )] d0 . (13) 



fc=0 



This is the well-known exact solution of the Ising model as found, e.g., in [17, 
24]. 



3.2 The thermodynamic hmit 

In order to facilitate the evaluation of /„ and we consider the limit 
oo, i.e. we replace 



1 



N 



k=0 



271 



2tt 



h{(l))d(f) . 



We label the functions obtained thereby by and g^, respectively. The 
error induced by taking the limit can be bounded using the following theorem. 

Theorem 2 (cf. [25]) Let h : W W be analytic and 27r -periodic. Then 
there exists a strip D = M x (—a, a) C C with a > such that h can be 
extended to a bounded holomorphic and 271 -periodic function h : D ^ C Let 
M = sup^|/i|. Then, 



e :- 



1 N~l ^ „27r 



/i(0)d0 



< 



2M 



It is straightforward to see that for any a > 0, the addends in both ( 1121) and 
(JT3ll are holomorphic on D = R x (—a, a) and bounded by 



M 



el"l"(l + e")(3 + exp[t(e 
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a/2 



-a/2 



)])• 



Applying the above theorem for a = 2.9, and taking into account that n < 
N/2, t > A, and N > 20, one obtains that |/n - /^| < e, - < e, 
where 

e < 30e-i-^^^+^-^* . (14) 
Taking the hmit in (fT2ll and (fT3l) . we find 



and 



/~ = — / [sin((2n - l)a) + sin((2n + l)a)] sin(2t sin a) da 
2vrJ 



[J2n-l(2t) + J2n+l(2t)] 





1 

2 

2nJ2„(2t) 



2t 



(15) 



1 

= — / [cos(2na) + cos((2r2 + 2)a)] cos(2t sin a) da + In 
27r Jo 

= ^[J2n{2t) + J2n+2{2t)]+In 

= + (16) 

where /q = |, /_i = — |, and /„ = otherwise. Here, Jk are Bessel func- 
tions of the first kind, and we have used the recurrence relation ^[Jn-i{z) + 

Jn+l{z)] = nJn{z)/z. 



3.3 Lower bound for the block entropy 

The reduced state of the first L sites is again Gaussian and its correlation 
matrix is given by the 2L x 2L upper diagonal subblock of Tt, which we label 
by At. Its entropy can be computed from the normal mode decomposition [26] 

L 

1=1 

1+x 1+x 1—x 1—x 

log log . 

2^2 2^2 



where O G S0(2L), as 



with 

h{x) = 



8 



We lower bound h{x) > 1 — a;^ as in [23,27] and thus obtain the bound 



SL{t) > E 1 - 
1=1 

= L + ltr[A',] . 

Writing 

A C 
-C^ B 

in the 1, . . . , L and L + 1, . . . , partitioning (thus A = At etc.) and using 
that the purity of the overall state implies Tf = —1 [22,23], we find A"^ — 
CC^ = -1 and thus 

We further bound the sum by only considering the lower left and the upper 
right corner of C, i.e., all entries 7„ in ffTUl) with l<n<L or N — L<n< 
N — 1 (equivalently, —L < n < —1, since ■jN-n = 1-n), and obtain 

SL{t)>lY. Hll7n||^ (IT) 



2 

n=—L 



We have 



=:A.n ='.Bji 

and with < l/-\/2 for \n\ > 1, 

S„<e(^4v^M + 2|/„|„||j , |n|>l. 
Thus, the total error made in ffT7|) is bounded by 

ltNB,<e(?^M!£±i).l)<0.3 
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where the latter follows using ffT^ together with L < N/2, 4 < t < and 
> 20. Thus, 



1 ^ 

Slit) >^Y. 1^1^" - 



- 2 

n=—L 



^ifjim Ji(2t) 1 



fc=l 



(2t)2 2t 4 

Pj|(2t) 



> V^l^^^ - 0.14 



A;=l 



where we have used Eqs. ^ and ([16]), J_„(^) = (-1)"J„(^), |Ji(2t)| < 
1/a/2, and t > 4. We now bound the sum in (1181) using Lemma [T] and [2] (see 
Appendix ) and obtain 

4 

> — t - ilnt- 1 

using 4 < t < eL/4 and L > 10, which proves Theorem [T] □ 
Note that different constraints on t, L, and can be chosen, resulting in 
different corrections to the leading terms in the expansion of Siit). 



4 Discussion 

We gave a rigorous proof to the observation that the entanglement entropy 
in a simple one- dimensional system increases essentially linearly. Although 
the proof is tailored to a single exactly solvable model, we think that this 
is a widespread property, thus being a requirement which has to be met 
by every classical simulation method. Of course, there always exist specific 
cases with tailored representations (e.g. the quasi-free Fermions used here or 
Clifford circuits) close to which one can meet this requirement [28,29] - even 
MPS allow for a description of specific states with only polynomial effort 
in the entanglement entropy (e.g., if the matrices are tensor products [30]). 
Yet, these methods seem much more restricted in their applicability than 
MPS, so that based on present knowledge it seems unlikely that there is 
an efficient and generally applicable classical simulation method. Note that 
similar results can be found by considering any Renyi entropy with a > 1 
instead of the von Neumann entropy [19]. Also observe that in our argument 
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we required a good global approximation of the state which is a priori not 
necessary for obtaining good results for observables with only local non-trivial 
support. 
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Appendix: Bounds on Bessel sums ^n^Jni^) 



Lemma 1 For z > 1, 



Q{z) ■.= Y,n^Jl{z)>^z'-\zHn 
^ — ' ovr 



n.=l 



4 - TT , Stt - 4 
z ■ — 2 



47r 



127r 



Proof. — Differentiating Q using the recurrence relations 

2J'^{Z) = Jn^l{z) - Jn+l{z) 



yields 



Jn{z) = —[Jn-l{z) + Jn+l{z)] 



n=l 



J2(z)+X^4n42(z) 



n=l 



To obtain a lower bound on f{z), we differentiate again and find 

CXj 

f{z) = 2zJ2[Jn-l{^) - Jlli^)] 
n=l 

= 2z[j',{z) + J',{z)] . 



As we prove in Lemma [3] later on. 



nz z 
11 



2 ' 



Z> 1 . 



(19) 



(20) 



(21) 



(22) 



Bounding Jq{z) + Ji{z) > for < z < 1, we obtain by integration 

/(;,)>liillll_21nz 

71 

for z > 1 and f{z) > otherwise. Together with J^iz) > 0, we have 
Q'{z) > zf{z)/2 and by integration, ( |T9l) follows for z >1. □ 

Note that the proof can be extended straightforwardly to sums J2'^=k- 



Lemma 2 For z < eK/A, K > 2, 

Kz"^ ^ ez\^K 

n-^j;[z) < 

n=K+l 



Qiz) = Y: < ^ (^) • (23) 



Proof. — As in the proof on Lemma [T], we differentiate Q using the recurrence 
relations (!20l) and find 

oo 

Q\z) = '-[{K+lfJl{z) + {K + 2fJl^,{z)+ J2 4r^^nW 

n=K+2 

' V ' 

and 

f'iz) = 2z[Jl^,{z) + Jl^,{z)] . 

We now use the lower bound |Jn(^)| < {z/2)"'/n\, which after integration 
leads to 

Qiz)<CiK,z) (-) 

with a prefactor C{K, z) of order {1/K\y, which for z < eK/A, K > 2, and 
using k\ > J^°° e~*t^dt > (k/e)^, straighforwardly reduces to (1231) . □ 



2K+2 



Lemma 3 For z > 1, 



Jl{z) + Jl{z)> 



TTZ Z^ 



Similar bounds can be proven for any pair + J^_^_i along the same lines. 
Note that the asymptotic scaling J^iz) + J^^i{z) ~ 2/7r2; is well known [31]. 
Proof. — We use the expansion ( [31], Sec. 7.3) 



Jn(,z) 




2 

TTZ 



cosU 



niT n , 



(n, 2m) 



\m=0 



{2z] 



2m 



P(z) 
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sm z 



\»n=0 ^ ^ 



with 

r(n + m + 1/2) 



m!r(n — m + 1/2) ' 

where the remainders P and Q he between zero and the first neglected term 
of the respective series, given that 2>0, 2p>ri — |, and 2g > n — |. With 
JniA — Sniz) + En{z), whcrc Sn{z) comes from the series and En{z) from 
the remainders, we have | Jn(2:)| > |5'n(^)| — |-E'n(-2)| and thus 

4{z)>Sn{zY-2\Sn{z)\\E^{z)\ . (24) 

We choose p = q = 1, and first upper bound |5'„(2;)| < (n = 0, 1) for 

3 



z > ^ and then, using ([2 

2. . 2/ N 2 1 3 1 45 2 1 

Jo (^) + ^1 (^) > — - ^ - ^^y^, + - 32^^^4 ^ ^ - 

for z >1, using | sin(2)| < 1, | cos(2;)| < 1. □ 
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